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ABSTRACT 

We present the data reduction pipeline for the Hi-GAL survey. Hi-GAL is a key project of the 
Herschel satellite which is mapping the inner part of the Galactic plane (|/| < 70° and \b\ < 1°), 
using 2 PACS and 3 SPIRE frequency bands, from 70/mi to 500/mi. Our pipeline relies only 
partially on the Herschel Interactive Standard Environment (HIPE) and features several newly 
developed routines to perform data reduction, including accurate data culling, noise estimation 
and minimum variance map-making, the latter performed with the ROMAGAL algorithm, a 
deep modification of the ROMA code already tested on cosmological surveys. We discuss in 
depth the properties of the Hi-GAL Science Demonstration Phase (SDP) data. 

Key words: instrumentation - data reduction 



1 INTRODUCTION 

The Herschel Space Observatory was launched from Kourou in 
May 2009 aboard an Ariane 5 rocket. Two of three scientific in- 
struments on the focal plane (PACS and SPIRE) are capable to ob- 
serve the infrared sky with unprecedented angular resolution and 
sensitivity, providing photometric observations in 6 d ifferent bands 
(70mn, 100/zm, 160//m, 250//m, 350//m and 500iUm: [Pilbratt et all 
l20ld and reference therein). 

The PACS photometer is composed of two bolometer arrays: a 
64x32 pixel matrix arranged from 8 monolithic subarrays of 16 x 16 
pixels each centered on the 70yt/m and 100//m wavelength (blue and 
green bands), and a 32 x 16 pixel matrix organize d in two subar- 
rays fo r the band centered on 160//m (red band), see Pog litsch et alj 

boioh . 

SPIRE comprises a three band photometer, operating in spec- 
tral bands centered on 250yt/m, 350yt/m and 500//m. Every band uses 
a matrix of germanium bolometers (139, 88 and 43 respectively) 
coupl ed to hexagonally packed conical feed horns dGriffin et al.l 
Hoifl). 

In order to handle science data provided by the Herschel in- 
struments, including the data retrieval from the Herschel Science 
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Archive, the data reduction through the standard pipelines and the 
scientific analysis, an official software environ ment called HI PE 
(Herschel Interactive Processing Environment, lOttetalJEoioh is 
available from ESA. 

The raw data provided by the satellite are reduced in HIPE to 
generate scientific data (so-called Level 2) and intermediate prod- 
ucts of the data reduction process (Level and Level 1 data). 

In this paper, we describe the dedicated pipeline created to 
obtain maps for Hi-G AL (Herschel Infrared Galactic Plane Survey, 
Moli nari et al.l2010ah . Hi-GAL aims to homogeneously cover with 
observations in 5 contiguous IR bands between 70yt/m and 500yt/m 
a 2 degrees wide stripe of galactic plane between / = -70° and 
/ = 70°. 

The Galactic plane shows emission varying from point-like 
sources to large-scale structures and with intensity varying over a 
wide dynamic range. In this work we show that the existing stan- 
dard reduction strategy (which is based on HIPE version 4.4.0, re- 
leased in November, 1 1th 2010) is not optimized to reduce Hi-GAL 
data and that a dedicated pipeline can enhance the quality of the 
Level 2 products. 

After Herschel successfully passed the Performance Verifica- 
tion Phase (PV Phase), two fields of the Hi-GAL survey were ac- 
quired during the Science Demonstration Phase (SDP): 2x2 square 
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degree areas of the Galactic plane centered on 30° of longitude 
(hereafter, / = 30°) and on 59° (/ = 59°). 

We describe the data reduction tools used to obtain high qual- 
ity maps from SDP data, with the aim to provide a reliable en- 
vironment for the Routine Phase (RP) data. The maps provided 
by our pipeline are suc c essfully used for several works like, e.g . 
. iMolinari et"ai1 (l2010bh . [Martin et all (l201Qh . |Peretto et al.1 feoioh . 

The paper is organized as follows: in Section [2] we describe 
the acquisition strategy for Hi-GAL data; in Section [3] we describe 
the pre-processing steps of the data reduction pipeline, necessary to 
prepare data for the map making and the tools that we have devel- 
oped to that purpose. In Section[4]we describe the ROMAGAL map 
making algorithm used in order to estimate the maps. ROMAGAL 
is used in place of the MadMap code which is the map making 
algorithm offered in HIPE. The quality of the ROMAGAL maps 
for both PACS and SPIRE instruments, related to the SDP obser- 
vations, will be analyzed in Section [5j in Section [6l we draw our 
conclusions. 



2 THE Hi-GAL ACQUISITION STRATEGY 

Hi-GAL data are acquired in PACS/SPIRE Parallel modeQ, in 
which the same sky region is observed by moving the satellite at 
a constant speed of 6(y7sec and acquiring images simultaneously 
in five photometric bands: 70/um and 160//m for PACS and 250//m, 
350yt/m and 500//m for SPIRE. 

The whole data acquisition is subdivided in 2° x2° fields of sky 
centered on the Galactic plane. Every Hi-GAL field is composed of 
the superposition of two orthogonal AOR (Astronomical Observa- 
tion Requests). Each of them is based on a series of consecutive, 
parallel and partly overlapped scan legs covering 2° x 2° square de- 
grees /Thescamn^gstotegy adopted for Hi-GAL is fully described 
in (Molin arietaflhoiOah . The superposition is performed in or- 
der to obtain a suitable data redundancy and for better sampling 
the instrumental effec t like the high-frequency detector response 
dMolinari et al.ll2010ah . 

The acquisition rate for the parallel mode is 40 Hz for PACS 
and 10 Hz for SPIRE, although the PACS data are averaged on- 
board for an effective rate of 5 Hz and 10 Hz for the lOfim and 
160yt/m array respectively. The implications of the PACS data com- 
pression are detailed in Section [53] 

An example of the scanning strategy of the Hi-GAL survey 
is shown in Figure Q] The coverage map of the PACS blue array 
is shown on the left panel and on the right panel we highlight the 
superposition of one scan leg to the following ones, by enlarging 
the bottom right corner of the left image. Two calibration blocks 
for each AOR, during which the instrument observes two internal 
calibration sources located on the focal plane, were provided during 
the SDP observations. They appear as higher than mean coverage 
areas and are marked by black and green circles for the 2 AORs 
in Figure \T\ Higher coverage zones are also clearly visible in the 
slewing region at the end of each scan leg, where the telescope 
decelerates and then accelerates before initiating the next scan leg. 

3 MAP MAKING PRE-PROCESSING 

The aim of the pre-processing is to prepare Hi-GAL data for map- 
making. 

1 http://herschel.esac.esa.int/Docs/PMODE/html/parallel_om.html 




Figure 1. Left: the coverage map of PACS blue array. Right: a zoom of 
the bottom right corner where is clear the effect of the superposition from 
one scan-leg to the next. Black and green circles highlight the calibration 
blocks that, during PV and SDP phase, were observed twice during which 
the internal calibration sources are observed. 



While map making is performed using a Fortran parallel code 
borrowed from cosmological observations (see Section |4), the pre- 
processing is done through a series of IDL and jython tools to be 
run on the data one after the other. After having tried to map Hi- 
GAL data using the standard tools provided within HIPE, we de- 
cided to develop our own routines that we tailored specifically for 
the reduction of data affected by bright and irregular background, 
as in the Galactic plane. In fact, high-pass filtering used in HIPE 
to cure the long-time drift also removes a large part of the dif- 
fuse Galactic emission. Furthermore, the standard deglitching em- 
bedded in HIPE , the MMT (Multiresolution Median Transform, 
IStarck et al .11 19951) . generates false detections in correspondence to 
the very bright sources when we apply this task to the PACS data. 
On the other hand, the deglitching procedure based on wavelet anal- 
ysis used by SPIRE does not affect the data, given also the lower 
spatial resolution compared to the PACS one. We therefore use the 
HIPE task for SPIRE data only. 

Herschel data are stored in subsequent snapshots of the sky 
acquired by the entire detector array, called frames. In a first pro- 
cessing step, the standard modules within HIPE are used to gener- 
ate Level 1 data for both PACS (except the deglitching task) and 
SPIRE. Thus, the data are rearranged into one time series per de- 
tector pixel, called Time Ordered Data (TOD), in which the cal- 
ibrated flux ( Jy beam" 1 for SPIRE and Jy sr 1 for PACS) and the 
celestial coordinates are included. At the end of this step, TODs are 
exported outside HIPE in fits format. In the subsequent processing 
steps, TODs are managed by a series of IDL tools, in order to pro- 
duce final TODs free of systematic effects due to electronics and of 
glitches and corrupted chunks of data due to cosmic rays. To each 
TOD a flag file is attached to keep track of any flagging done during 
data reduction steps. 

Preprocessing includes identification of corrupted TODs (or 
part of them), drift removal and deglitching. The following steps 
will be the Noise Constraint Realization (NCR) and ROMAGAL 
mapmaking; they will be both described in detail in the next Sec- 
tions. The summary of the entire pipeline is shown in the diagram 
in Figure |2] 

3.1 Corrupted TODs 

A TOD can be partially corrupted by the random hiting of charged 
particles (cosmic rays) which produce strong and spiky signal vari- 
ations called glitches. Two different effects on the TODs can be 
identified: the glitch corrupts a single or few consecutive samples, 
generating spiky changes along the TOD. This is the most common 
effect and in Section [331 we describe how to detect and mask the 
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Figure 2. Schematic representation of the Hi-GAL data reduction pipeline 

data for PACS, as well as to mask any possible residual glitches 
for SPIRE. Very powerful glitches, on the other hand, make the de- 
tector signal unst able for a considerable interval of time (see, e.g., 
Billot et alJ2010h . This effect depends on the glitch amplitude and 
on the bolometers time response. These events affect a much larger 
part of the TOD that cannot be used for mapmaking. 

Their impact results in a bias on the bolometer timeline with 
a low-frequency drift which can involve a considerable number of 
samples, as shown in Figure [3] 

In that Figure, blue crosses represent the observed timeline of 
one bolometer of the blue array. 

Automatic identification of the (partially) corrupted TODs 
exploits the first derivative of the TOD to detect extraordinary 
"jumps" in the signal. In order to determine what portion of the 
TOD is to be flagged, the flagging tool exploits the fact that the de- 
tector pixels response that have been hit by a cosmic ray is mostly 
exponential. Data samples ranging from the jump to the time sam- 
ple at which an exponential fit reaches 98% of the signal level be- 
fore the event are identified as bad data and stored as such in the 
corresponding flag file. In case the exponential fit does not reach 
98% of the starting value before the end of the TOD, then all data 
starting from the hit will be flagged as is the case in Figure |3] This 
procedure is applied both in the cases of a changing in responsivity 
or a detector offset alteration. In the latter, we estimate the fit with 
an exponent equal to 0. The described procedure was adopted to 
process both PACS and SPIRE data. 

3.2 Drift removal 

After having identified corrupted data we proceed to the elimina- 
tion of changes in responsivity over time. The procedures are in 
principle identical for PACS and SPIRE data, the only differences 
account for different composition of the detector arrays and the data 
acquisition of the two instruments. 

The signal in PACS TODs that exit HIPE does not represent 
the true sky but is dominated by the telescope background and the 
(unknown) zero level of the electronics. The electronics further in- 
troduce significant pixel-to-pixel offsets. For each TOD, we miti- 
gate the effect of pixel-to-pixel offset by calculating and then sub- 
tracting the median level for each pixel from each readout. This 
ensures that all pixels have median level equal to 0. The median 
value is preferred over mean because the median is a much better 




Frames 



Figure 3. Timeline of a PACS blue bolometer. The exponential decay illus- 
trates the change in responsivity after frame 40000 due to the impact of a 
powerful glitch. 




Figure 4. Blue PACS frame after the median subtraction on each pixel. 
Diffuse emission and compact source are visible in the frame. 



representation of the sky+telescope background flux, and is much 
less sensitive to signal from astrophysical sources. 

Also, the subtraction of this offset from each TOD does not al- 
ter the signal in the final map, but it introduces only a global offset 
constant over the entire area covered in the observation. However 
it should be kept in mind that bolometers are inherently differential 
detectors which bear no knowledge of an absolute signal value; be- 
sides, any optimized map making methods like the one we employ 
(see Section [4} produce maps with an unknown offset value which 
needs to be independently calibrated. So it is important to reduce all 
the bolometers to the same median value, regardless of its amount. 
The pixel-to-pixel median subtraction has the effect seen in Fig- 
ure!?] Diffuse emission and compact sources are clearly visible in 
the frame. 

Still, when plotting a detector pixel timeline we see that the 
signal decreases systematically from the start to the end of the ob- 
servation. This trend is likely due to a combination of small changes 
in the thermal bath of the array and to small drifts in the electronics. 
The former affects the entire array, while the latter affects subunits 
of the detector array (in fact, PACS blue is divid ed into 8 units, 
PACS red into 2, electronically independent units (Po ditsch et al.l 
2008)). The drift is then a combination of these two effects: drifts of 
the entire array and drift of a single subunit. These effects are dom- 
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Figure 5. Median behavior computed on the whole array for each frame. 
The slow-drift behavior is due to the electronics and the thermal bath. 



inant with respect to the 1 // noise pattern, which will be described 
in the next Sections. 

It is in principle not a trivial task to decide which drift has to be 
subtracted first: the drift from the thermal bath (affecting the entire 
detector array) or the drift from the readout electronics (affecting 
sub-arrays differently)? Ideally both should be subtracted, if only 
it were possible to separate each component accurately, as the net 
effect in the data is the sum of both. 

Our methodology for removing the correlated signal drifts (on 
both the bolometer module/unit level and the array level) is based 
on tracing the low signal envelope of the unit or array median lev- 
els. In Figure [5] this envelope is the curve defined by the lowest 
signal values, and estimated as follows: 

i We compute the median value of the entire bolometer ar- 
ray/unit for each bolometer readout. Figure [5] shows one example 
for the entire array. 

ii The median values thus obtained are segmented and grouped 
by scan legs. Each scan leg is composed of ~ 1000 frames and we 
observed 54 scan leg for each 2° x 2° Hi-GAL field. 

iii For each scan leg we compute the minimum value of the ar- 
ray/unit medians. 

iv The resulting set of minimum median values for all scan legs 
are fit with a polynomial. 

The median value for each array/unit readout is chosen be- 
cause it is closest to the actual sky+telescope background. How- 
ever, as clearly seen in Figure [5] in the presence of strong astro- 
physical sources the median value is incorrect for our purposes. 
The strong sources appear as signal spikes in Figure [5] Hence, we 
take the next step of finding the minimum value from the set of me- 
dians belonging to a single scan leg. The idea is that at some point 
during the scan leg the median was indeed a true representation of 
the local sky+telescope and relatively free of source emission. This 
step allows us to reject the sources at the expense of degrading our 
time-resolution to scan-leg duration (~ 240 sec). The polynomial 
fit allows us to estimate the drift behavior at the same time resolu- 
tion as the PACS signal (5Hz and 10Hz for 70/um andl60//m band 
respectively). We further note that the correlated signal drift is rel- 
atively flat over a single scan leg; hence, the minimum value is not 
significantly affected by the presence of the monotonic signal drift 
itself in the first place. 

The minimum median method discussed above removes back- 
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Figure 6. Interpolation of the minima of the median evaluated on every 
scan leg. Each curve refers to a PACS subarray. The curves mimic the same 
behavior with slopes due to the different subarray electronics. 
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Figure 7. Minima of the median evaluated on the whole PACS blue array 
after the subarray drift subtraction. The dispersion is due to the intrinsic 
efficiency of the drifts subtraction tool. There is no residual behavior and 
the scatter is one order of magnitude under the value of the initial drift. 



ground levels from spatial emission structures that are of the or- 
der or larger than the scan legs yet preserves the spatial structures 
smaller than scan legs . In essence, information about the a bsolute 
calibration zero-point dMiville-Deschenes & Lagachell2005l) is lost 
but all spatial structures within our map boundaries are preserved. 

In Figure [6] we reported the minimum median values of each 
subarray. The common downward trend is due to the (common) 
thermal bath drift, while the different readout electronics are re- 
sponsible for the differences in the subarray slopes. 

We therefore decide to subtract the subarray drifts in order to 
consider both the thermal bath and the readout electronics behav- 
iors, but separately for each subunit. 

Once the individual subarray drift is removed, the remaining 
dispersion on the whole array is only a residual scatter due to the 
intrinsic efficiency of the removal tool, as shown in Figure |7] 

SPIRE array detectors are not divided into subarrays, so every 
procedure that has to be run 8 or 2 times for PACS data is only per- 
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Figure 8. Blue curves: interpolation of the minima of the median for one 
subarray of PACS blue channel when the scan direction is almost parallel 
to the Galactic plane. Red line: the fit that we choose to evaluate the drift 
without considering the minima affected by the Galactic emission (central 
bump). 



formed once per SPIRE band. SPIRE uses blind bolometers (5 in 
total) as thermistors to evaluate the most relevant correlated noise 
component: the bath temperature fluctuations. A standard pipeline 
module uses this information to perform an effective removal of 
the common drift present along the scan observation. HIPE also 
corrects for the delay between the detector data and the telescope 
position along the scan, using an electrical low-pass filter response 
correction. But despite these initial (and very effective) corrections, 
we apply the drift removal tool to SPIRE data in the same way as 
for PACS data: we fit a polynomial to the minimum of the me- 
dian of each scan leg (calculated over the entire detector array), 
that we then subtract to all TODs. Experience on the data shows 
that a residual long time drift is often present in SPIRE data. 

Finally, when removing drifts it is important to know how the 
observational scans are oriented. In fact, as the Galactic plane is 
very bright, scans across the plane will give rise to an increase of 
signal, on top of the general drift. On the other hand, when the scan 
is almost parallel to the plane of the Galaxy, the signal can domi- 
nated by its bright emission, also on the evaluation of the minima 
of the median. 

In this case, the curve to fit is estimated only in the scan legs 
where the median is not affected by the signal. 

Since the procedure is not automatic, care has to used when 
choosing what polynomial to fit and subtract from the data, in order 
not to eliminate genuine signal from the sky. From our experience 
the best choice is a first or a second degree polynomial, depending 
on the signal behavior observed. 

Higher polynomial degree can be necessary when part of the 
drift has to be extrapolated in order to avoid signal contamination. 

An example for one subarray is shown in Figure [8] 

Each scan leg is bigger then the overlapping region of the two 
scan directions for the scanning strategy adopted by Hi-GAL (see 
FigureQ])- Since the Hi-GAL fields are squared regions, the slowly- 
traversed direction of the AOR within the overlapping region have 
a length comparable with the scan leg. Thus, we assume that even 
if there is a signal gradient along the slowly-traversed direction of 
the AOR, it is not filtered out by the array medians subtraction. 



3.3 Deglitching 

To remove outliers in the time series of the bolometers we exploit 
the spatial redundancy provided by the telescope movement which 
ensures that each sky pixel of the final map is observed with differ- 
ent bolometers. Outliers detection is done with the standard sigma- 
clipping algorithm: given a sample of N values, first estimates for 
the mean and for the standard deviations are derived; then all the 
values that differ from the mean by more than n standard deviations 
are considered outliers and removed. 

For this algorithm the choice of n, the parameter that defines 
the threshold above which a value is considered an outlier, is usu- 
ally arbitrary: a certain n is chosen, very often equa l to 3, without 
any statistical justification. Recently [Pezzutol J2010h has derived a 
formula that starting from the properties of the error function for 
a Gaussian distribution and exploiting the discreteness of the num- 
ber of available measures, relates n to the size of the sample. The 
formula is 

n = -0.569 + V-0.072 + 4.991og(A0 (1) 

As a consequence, in the central region of the map, where the 
coverage (and so n) is high, the number of standard deviation is 
larger than in the outskirt of the map where the coverage is low. 
For instance, if a sky pixel has been observed with 40 bolometers, 
the above formula gives n = 2.25; so, once we have estimated the 
mean m and the standard deviations cr, all the values x t such that 
ABS(x; - m) > 2.25<x are flagged as outliers. If a pixel has been 
observed with 20 bolometers the threshold lowers to 1.96<x. 

This procedure is automatically iterated until outliers are no 
longer found. However, the procedure converges within 1 iteration 
in ~98% of the cases in which we have applied the analysis. 

The outliers detection is done in this way for both instruments, 
however for SPIRE, as explained before, we also make use of the 
standard deglitching algorithm (wavelet based) implemented in the 
official pipeline. But we found some weak glitches left in the SPIRE 
TODs so that we decided to run our deglitching algorithm also on 
SPIRE data. 

The number of glitches found is on the average about 15%, a 
value which is likely larger than the real percentage. For PACS we 
are now working on a different way to associate each bolometer to 
the sky pixels, taking into account the finite size of the sky pixels. 
For the first test cases we run, the percentage of detected glitches is 
now around 5-6%. 



4 THE ROMAGAL MAP MAKING ALGORITHM 

The ROMAGAL alg orithm is b ased on a Generalized Least Square 
(GLS) approach (Lupton 1993). Since the TOD is a linear combi- 
nation of signa l and noise, we can model our dataset d k for each 
detector k as (Wright 1996): 

d k = Pm + n, (2) 

where P is the pointing matrix, which associates to every sample 
of the timeline a direction in the sky, m is our map estimator of the 
"true" sky and is the noise vector. 

The observed sky, Pm, is the map estimator of the "true sky" 
convolved with the instrumental transfer function and the optical 
beam. However, in case of circularly symmetric beam profile, m is 
a beam smeared, pixelised image of the sky. 

In this case the pointing matrix has only a non-zero entry 
per row corresponding to the sky pixel observed at a given time. 
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Since the beam profi les for PACS dPoditsch et al J2008h and SPIRE 
(Griffi n etalJl2008h are only weakly asymmetric we can put our- 
selves in this simple case. Note that the transpose of the P operator 
performs a data binning (without averaging) into the sky pixels. 

Equation [2 holds only if the noise vector of each detector 
is composed of statistical random noise, with Gaussian distribu- 
tion and null average. All the relevant systematic effects (offset, 
glitches) have then to be removed with an accurate data prepro- 
cessing before map production, as explained in Section [3] 

The formalism can be easily extended in the case of multide- 
tector analysis. In this case the vector d contains the data relative to 
each detector. Rather, one has to take care to upgrade also the noise 
vector n, accordingly with the correct noise value for each detector. 

The GLS algorithm produces minimum noise variance sky 
maps. Noise properties for each detector have to be previously es- 
timated and provided in input to the algorithm as described in Sec- 
tiongj] 

The GLS estimate for the sky, m, is dNatoli et al.ll200ll) 

m = (P r N" 1 P)" 1 P r N" 1 d (3) 

where N = (nn r ) is the noise covariance matrix, which takes into 
account noise time correlation between different samples. Such 
correlation is particularly high at low frequencies because of the 
l/f (or long memory) noise. In case of uncorrelated noise (or 
white noise) the N matrix becomes diagonal and the problem is 
greatly simplified. If we further assume stationary uncorrelated 
noise, Equation [3]reduces to: 

m = (P T py 1 P T d. (4) 

P T P is the number of observations of a pixel of the map, so we are 
averaging the different TOD values into that pixels assigning the 
same weight to each sample. We will refer to this map estimate as 
"naive" or "binned" in the following. 

When non negligible noise correlation is present, as in the cas e 
of PACS JPoditsch etaT]l2008h and SPIRE dSchulz et alJ l2008h . 
Equation [3] must be solved. This is a challenging computational 
task since it requires, in principle, the inversion of the large (of the 
order of the number of pixels in the map) ma trix P r N -1 P, which 
is the covariance matrix of the GLS estimator (Lupton 1993). One 
key simplifying assumption is to invoke that the noise is stationary. 
In this case, the N matrix has a Toeplitz form which can be approx- 
imate ly treated as circulant, ignoring boundary effects dNatoli et al.1 
l200lt) . A circulant matrix is diagonal in Fourier space and its in- 
verse is also circulant, so the product between N _1 and a vector is 
a convolution between the same vector and a filter provided by any 
of the rows of the matrix. In the following we will refer to any of 
these rows as a noise filter. Its Fourier transform is the inverse of 
the noise frequency power spectrum. 

Considering the conditions listed above, the GLS map making 
algorithm performs the following operations, starting with rewrit- 
ing the Equation [3] in the form 

CP r N _1 P)mo - P r N _1 d = r (5) 

where m is the starting map used at the first iteration, gener- 
ally the naive map. 

The Pm product projects the map onto a timeline. Applica- 
tion of N _1 : this is a convolution which can be performed in Fourier 
space. Application of P T : this step projects the convolved timeline 
back into a map. 

The second term performs the convolution with the filter (ap- 
plying N _1 to the data vector d in Fourier space) and then the pro- 



jection of the convolved timeline into a map (applying P T to the 
product N _1 d). 

Then, we need to evaluate the residual r. If the residual is 
higher than a fixed threshold, it is used to p roduce a new map, 
mi , as described in Hestenes & Stiefell dl952h . This map will be 
considered instead of m for evaluating again the Equation [5] until 
convergence. This is achieved by running a Conjugate Gradient al- 
gorithm, an iterative method useful t o obtain a numerical solution 
of a system dHestenes & Stiefell 1952h . until convergence is reached 
with the residual lower then the threshold. 

The algori thm outlined is the same described in ("unroll, con- 
volve and bin": iNatoli et al.ll200ll) and is implemented in the RO- 
MAGAL code. The next section explains the strategy employed 
to estimate the noise filters used by ROMAGAL directly from in- 
flight data. 

4.1 Noise estimation 

In order to estimate the noise filters for ROMAGAL we need to 
investigate the noise statistical properties in the timelines. Data 
are mostly affected by two kind of statistical noise: l/f noise due 
both to the electronics and thermal background rad iation from the 
telescope or the instrume nts, and photon noise fsee lPoglitsch et al.1 
l2008l : ISchulz eUd1l2008h . 

The detector l/f noise arises in the electronic chain and it 
impacts particularly regions with low signal-to-noise ratio (SNR), 
where only diffuse emission is present. In those regions it can be of 
the same order of magnitude of the signal or even higher. In these 
cases GLS treatment is particularly effective. 

Photon noise is due to statistical fluctuation in the photon 
counts. This process follows poissonian statistic, so the SNR is 
proportional to the square root of counts. Since poisson distribu- 
tion tends to be Gaussian for large numbers, we can approximate 
photon noise as Gaussian on the map if the number of counts is 
large enough. 

Since bolometers are organized in matrix and sub-matrix, the 
signal of a bolometer can be correlated with the signal of another, 
generally adjacent, bolometer. These effects could be both statis- 
tical and deterministic. We already described how to remove the 
deterministic common mode from TOD (like the thermal bath vari- 
ations, see Section [3}. 

One possible source of statistical cross-correlated noise is the 
crosstalk between bolometers: the signal of one pixel may contam- 
inate the signal of its neighbors through the capacitive inductive 
coupling, generating a common mode called "electrical crosstalk". 
On the contrary, "optical crosstalk" is due to diffraction or aber- 
rations in the optical system which co uld drive an astron omical 
source to fall on inappropriate detectors (IGriffin et al . 2008). 

We then analyze the residual contribution of the statistical 
component of the correlated noise. We found that the residual cor- 
related noise level in each pixel is negligible with respect to the in- 
trinsic detector noise level for both PACS and SPIRE instruments, 
as described in the following. 

In principle, the noise properties vary significantly across the 
array and we had to estimate the noise power spectrum for each 
bolometer. To do that we have processed "blank" sky mode (i.e. 
filled with negligible contribution from sky signal) data acquired 
during the PV Phase. 

In Figure [9] we show a typical noise spectrum estimated for a 
pixel of the 160//m PACS band (black) and the cross-spectrum be- 
tween two adjacent bolometers (red). The cross-spectrum evaluates 
the impact of the cross-correlated noise in the frequency domain 
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Figure 9. Black line: typical noise spectrum of a PACS 160//m detector, 
estimated on blank sky data. Red line: cross spectrum between two detectors 
of the same subarray. The level of the cross-correlated noise is significantly 
under the noise level of each single bolometer, so we can reasonably neglect 
it. 
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Figure 10. Same as Figure |9] for a SPIRE 250//m bolometer. For SPIRE 
also the noise level of cross-spectrum is reasonably negligible with respect 
to the auto spectrum level. 

between two different bolometers. The level of the cross-correlated 
noise is at least 4 order of magnitude below the auto-correlated 
noise power spectrum of each pixel. Note that this means we do 
not see any relevant cross-correlated noise, despite the fact that 
crosstalk can be present into the timeline. 

In Figure [TO] we show noise power spectra of the 250//m 
SPIRE band bolometers. Also in this case the cross-spectrum is 
negligible. 

Noise spectra of both PACS and SPIRE display low-frequency 
noise excess (1//). In case of SPIRE spectra (Figure fTOl a high fre- 
quency rise is also evident, which is due to the deconvolution of 
the bolometric response function. PACS spectra do not show this 
behavior because the bolometer transfer function is not deconvo- 
lved by the standard pipeline. 

4.2 From ROMA to ROMAGAL 

The ROMAGAL GLS code has been optimized to recover the fea- 
tures in the Hi-GAL fields with high accuracy. 

Hi-GAL observes the Galactic plane where the dynamic range 



Band 


Total Time (sec) 


RAM (Gb) 


70^ 


~ 1400 


16 


160//m 


~ 1000 


8 


250//m 


~ 180 


4 


350//m 


~ 130 


1 


500//m 


~ 100 


1 



Table 1. Time and minimum RAM amount required from ROMAGAL for 
each PACS and SPIRE band using 8 BLADE processors. 

of sources spans over several orders of magnitudes. This poses 
strong constraints on the map making algorithm: both the weak dif- 
fuse emission and the bright compact sources in e.g., star forming 
regions have to be recovered with high accuracy. The signal often 
exhibits steep gradients that are hard to follow for the GLS solver, 
which relies on the assumption that the sky signal does not vary sig- 
nificantly within a resolution element (see below l53l >. At the same 
time, several systematics affect the dataset. As explained above, 
many of them are cured at the preprocessing level. However, their 
removal generates a conspicuous amount of transient flagging, that 
must be correctly handled by the GLS code. 

The core of ROM AGAL is the same of the ROMA code 
jde Gasperis et al.l2005l) where the input-output routines have been 
deeply modified to adapt to the HIPE generated dataset. ROMA- 
GAL inputs are the TOD generated by HIPE, pointing and transient 
flagging. These have the same format for both PACS and SPIRE. 
ROMAGAL outputs are fits file containing the optimal map, ar- 
ranged in standard gnomonic projection routines. The code is writ- 
ten in FORTRAN 95 and relies on the MPI library for parallel calls. 
It runs on the Hi-GAL dedicated machine, called BLADE, a cluster 
of 104 processors at 2.5 GHz each and 208 Gb RAM total. Its nodes 
are interconnected with MPI-infiniBAND. The machine is located 
at IFSI-Rome. 

As explained in the previous Section, the computation of 
Equation [3] is unfeasible due to the size of the map's covariance 
matrix. However, we assume the noise of each Hi-GAL field to be 
stationary to set up an FFT based solver built upon a conjugate 
gradient iterative algorithm (see Section |4). Such a scheme can es- 
timate the final maps with a precision of order of e = 10~ 8 in ~ 150 
iterations for Hi-GAL . ROMAGAL computational time scales lin- 
early with the size of the dataset and only weakly with the number 
of pixels in the output maps. The scaling with the number of proces- 
sors is highly optimal in the range of cores required for the Hi-GAL 
analysis (< 50). For the largest channels (PACS blue band), a final 
GLS map of a 2° x 2° field requires about 16 Gbytes of RAM and 
1400 sec on 8 cores. Due to the high number of array pixels (2048), 
this channel is the largest dataset to analyze as well as the most de- 
manding in terms of computational resources. Further information 
on resource consumptions can be found in Table Q] 

4.3 Optimal treatment of transient flagging 

As mentioned above, the timelines are inspected for bad data sam- 
ples that must be excluded from map making as part of the pre- 
processing pipeline. Bad data can arise due to a variety of reasons. 
They are generally caused by transient events, either unexpected 
(e.g., glitches, anomalous hardware performance) or expected (e.g., 
detectors saturating because of a bright source, observation of a cal- 
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ibrator). Once identified, a flag record is generated and stored for 
these anomalous events, so that their contribution can be safely ex- 
cluded from map making. Flags in the TOD pose a potential prob- 
lem for ROMAGAL because its solver is based on the FFT as dis- 
cussed in the previous section. The FFT requires the timeline to be 
continuous and uniformly sampled. Since noise in the PACS and 
SPIRE data is correlated, just excising the flagged sampled to fic- 
titiously create a continuous timeline would interfere with noise 
deconvolution, and is thus not a safe option. Instead, we advocate 
using a suitable gap-filling technique. The rest of this section is 
mostly devoted to defining which, among the various options for 
gap filling, is best suited for the Hi-GAL maps. 

It is important to realize that the output map will depend on the 
content of the flagged sections of the TOD even if these values are 
not associated to any (real) map pixel. This is due to the convolution 
performed within the solver that "drags" out data from the flagged 
section of the timeline, even if, when the data are summed back 
into a map, the P operator is applied only to the unflagged samples. 
Since one is not in control of the content of the flagged section of 
the timeline, a kind of gap-filling must be performed in any case. 
We have tested different recipes for gap-filling making extensive 
use of simulations. We have treated separately the signal and noise 
components of the timelines, running separately noise dominated 
and signal dominated cases, because the behavior towards flags of 
the two components is different, as it will be shown below. 

The simplest form of gap filling is to remove the content of 
the flags altogether, replacing the flagged sections with a nominal 
constant (usually null) value. This works very well on a signal- 
only simulation of the Hi-GAL field. However, it fails dramatically 
when noise is present, as evident from Figure [TT](left panel), where 
a markedly striped pattern in the reconstructed map is seen (in this 
simulation, the Hi-GAL noise has been amplified by a factor 100 to 
make its pattern more evident). The reason for this failure is read- 
ily understood: The GLS map making employed requires noise sta- 
tionarity (see S ection |4~2l above) , which is obviously not preserved 
by zeroing the gaps. A less obvious effect is that even if gaps are 
filled with a noise realization with the correct statistical properties, 
but unconstrained, the GLS map making is bound to fail as well, 
as shown in the middle panel of Figure [TT] A noise realization is 
said to be constrained wh en it respects certain boundary conditions 
(Hoffman & Ribak 1991), which in our case are represented by the 
noise behavior outside the gap. Unconstrained noise inside the gap, 
despite having the correct statistical properties, creates a border dis- 
continuity that cause s the GLS map maker to behave sub-optimally 
dStompor et al.ll2002b . We have employed a code to create Gaus- 
sian noise c onstrained (NCR) realiza tions, based on the algorithm 
set forth bv lHoffman & Ribald Jl99lh . The code uses in input the 
noise pattern and statistical properties, as measured from the time- 
lines. The results on noisy simulations are excellent, as set forth by 
the third (rightmost) panel in Figure [TT] Note, however, that Figure 
[TT]refers to a noise dominated simulation. We now turn to discuss 
the effect of a non negligible signal contribution (a far more realis- 
tic case for Hi-GAL ). 

We have verified that the presence of non negligible signal in 
the timelines does not affect the results provided that the NCR is 
performed using the underlying noise field as a baseline. Measur- 
ing the latter in presence of signal is however impractical. It would 
be significantly simpler if the NCR could be run directly on the 
timelines themselves, constraining thus the (fake) noise realization 
within the gap to the outside noise plus signal (true) values. Unfor- 
tunately, this poses a problem for Hi-GAL data: large signal jumps 
are present in the field and the resulting gap filling realizations are 




Figure 12. Top row: For a signal dominated case, shown are the relative 
differences between the input simulated map and the output obtained with 
ROMAGAL (in the "virtual map" mode) with NCR performed on the un- 
derlying noise (left), on the original signal plus noise TOD (middle), and 
without NCR after replacing the gaps with null values. The latter case is 
clearly striped, but no signal related artifacts are present. Bottom row: For 
a noise dominated case, shown again are the relative difference versus the 
input map obtained by ROMAGAL with "virtual map", assuming NCR on 
underlying noise only (left), without NCR (middle) and with NCR on the 
signal plus noise TOD (right). As expected, the left panel shows the best 
residuals, but the right one appears as a good compromise (see also text). 



affected in a non negligible manner by the boundary conditions, 
at least with the present version of ROMAGAL. This behavior is 
different from what happens for ex periments aimed at the Cosmic 
Microwave Background (see e.g. iMasi et al ] (l2006h ) where NCR 
codes are routinely run on the timelines as they are. In order to find 
a workaround that would spare us the inconvenience of estimating 
the underlying noise field to serve as a NCR input, we have modi- 
fied the flag treatment of the ROMAGAL code as explained in the 
following. 

The original version of ROMAGAL makes use of a single ex- 
tra pixel (dubbed "virtual pixel") to serve as junk bin where the 
contents of the gaps are sent when applying the P T operator within 
the ROMAGAL solver. This approach, as stated above, works ex- 
cellently in presence of both signal and noise, irrespective of their 
relative amplitude, provided the NCR code assumes the underlying 
noise field as a baseline to perform the realization. In order to relax 
this assumption, we have modified ROMAGAL to take into account 
not a single virtual pixel but an entire virtual map. In other words, 
we introduce a virtual companion for each pixel of the map, and use 
it as a junk bin to collect the output from the gaps they correspond 
to. The hope is to redistribute the content of the flagged sections 
more evenly, preventing artifacts. This approach obtains satisfac- 
tory results when the NCR code is run on the original (signal plus 
noise) timelines, as shown in Figure [12] 

To summarize our findings: 
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Figure 11. Shown are the results obtained in a noise dominated regime (normal Hi-GAL noise is amplified by a factor 100). Left panel is the map obtained 
replacing flagged data samples with null values (clearly it does not work), middle panel is the map obtained replacing data samples with unconstrained noise 
realization (does not work either), right panel shows the map obtained using our NCR code (does work). 



• Using a NCR realization code is always necessary in order to 
avoid artifact striping effect into the map 

• If the NCR code is run on the underlying noise-only timelines 
(which are cumbersome to estimate) we obtain the best quality out- 
put map, with no signal related artifacts and a noise standard devi- 
ation which is lower (by a factor ~ 2) with respect to the case in 
which no NCR is performed. 

• Running the NCR on the original timelines is possible with 
the "virtual map" approach. No signal artifacts are detected in the 
difference maps and the advantage in terms of noise standard devi- 
ation with respect to the no NCR case is still present, and of order 
10% to 20% on average. 

We have therefore chosen this latter approach as the baseline for 
the Hi-GAL pipeline. 



5 ROMAGAL MAPS 

In this section we analyze the final map obtained running our dedi- 
cated pipeline. We analyze the Point Spread Function (PSF) for the 
five bands in order to fix the resolution of the ROMAGAL maps. 
We compare the final GLS map with the naive map and we point 
out the differences and the capability to recover the diffuse emission 
from the GLS map. Finally we discuss about the noise residuals on 
the maps. 

5.1 Point Spread Function and pixel size 

The angular resolution (pixel size) of the final map is a free param- 
eter. Its choice is a compromise between competing requirements. 
A small pixel size assures a correct sampling of the PSF; 

indeed, assuming a Gaussian profile for the PSF (which is rea- 
sonable as discussed in the following), the Nyquist theorem im- 
poses that to better sample a 2-d image we need to set a pixel size 
which is at most one third of its FWHM value. 

On the other hand, a too small pixel size can cause the loss of 
redundancy, useful to reduce the white noise level, and (even) some 
non-observed pixels in the final map. 

The diffraction limited beam of PACS at 70/im is 5.2". Thus, 
we should build the map with a pixel size of at least 1.8". How- 
ever, due to the limited bandwidth available for transmission, es- 
pecially in PACS/SPIRE Parallel mode, PACS frames are coadded 



on-board the satellite before broadcasting. For the 10/im channel, 
a coaddition of 8 consecutive frames is applied by on-board soft- 
ware. Since the acquisition rate is 40Hz and the scanning speed for 
Hi-GAL is set to 60"/ s, two close frames are snapshots of the sky 
acquired 1.5" apart. Due to coaddition, the satellite provides one 
averaged frame every 12"; in spatial coordinates, this is twice the 
beam width of the PACS blue channel. The measured PSF then is 
not the diffraction limited beam but it results in an elongation along 
the scan direction due to the co nvolution of the coaddition with the 
beam. As shown in Eutzl 12009). the observations of Vesta and QfTau 
with the blue channel evidenced a FWHM equal to 5.86" x 12.16" 
as a result of a 2-d Gaussian fitting, elongated in the scan direction. 

PACS 160yL/m is also affected by the averaging on-board, but 
only 4 frames are coadded together. The nominal instru menta l 
beam is 12.0", while the measured is 11.64" x 15.65" jLutzl2009h . 
elongated along the scan direction. However, in this case we can 
sample the beam without issues, and the effect of coaddition on the 
final map is negligible. 

For the Hi-GAL fields the scanning strategy consists of two 
orthogonal AORs, therefore the redundancy regularizes the PSF, re- 
sulting in approximately 2-d symmetric Gaussian profile, as shown 
in Table |2] 

According to the values reported in Table |2] we observe quasi - 
symmetric beams with an averaged ellipticity of less than 15% for 
both blue and red channel and the axis oriented randomly with re- 
spect to the scan direction. 

We choose a pixel size of 3.2" for the PACS 70//m band 
which samples the observed beam at Nyquist frequency. Below this 
threshold, the diffuse emission areas become too noisy due to the 
low SNR. Similarly, we can choose a pixel size of 4.5" for red band 
without loosing redundancy. 

SPIRE does not suffer from on-board coadding and the de- 
tectors were built to reach the band diffraction limit. In-flight data 
show that the SPIRE beam is well approximated by a 2-d asym- 
metric Gaussian curve with the major axis orientation independent 
of the scan direction, with an ellipcticity not bigger then 10% (see 
Sibtho rpe et all (l2010h ). We set the pixel size for each SPIRE band 
equal to one third of the nominal beam. In Table |2] we also report 
the beam measured in the SPIRE maps. 

The average e llipticity we observe agrees with found by 
ISibthorpe etaPfcOlOh. On the contrary while the FWHM for the 
two axis found by Sibthorpe et al 1 (12010) are in agreement with the 
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Band 


Nominal Beam (arcsec) 


Measured Beam (arcsec) 


EUipticity 


Pixel Size (arcsec) 


70yt/m 


5.2x5.2 


~ 9.7x ~ 10.7 


14.6% 


3.2 


160//m 


12.0 x 12.0 


~ 13.2x ~ 13.9 


14.7% 


4.5 


250//m 


18x18 


~ 22.8x ~ 23.9 


8.3% 


6.0 


350//m 


24x24 


-29.3x31.3 


8.8% 


8.0 


500/im 


34.5 x 34.5 


~41.1x -43.8 


9.7% 


11.5 



Table 2. Nominal (2nd column), map-measured beam (two AOR, 3rd column) and ellipticity (4th column) of each band. 




Figure 13. Left: particular of the ROMAGAL map of the red PACS array, 
/ = 30° field. Right: particular of the naive map, same band and field. The 
l/f noise is evident on the naive map, as well as its minimization is evident 
on the GLS map without loosing the signal of the diffuse component. 



nominal (within the error), our measured beam results in a FWHM 
larger then the nominal of ~ 25%. 




Figure 14. Particular of map difference between PACS 160/mi igls and 
naive, same region of the previous Figure. The stripes due to 1 // removal 
made in the igls are evident. 



5.2 Hi-GAL SDP results 

The quality of the outcome can be qualified by the comparison be- 
tween the igls maps with the naive maps. In fact, the naive map is 
the simple averaging of signal recorded by every spatial pixel and 
it represents "the least compromised view of the sky". 

Since the TOD are created at the end of the preprocessing 
steps, when the data are a combination of only signal and 1 // noise, 
we expect the 1 // residuals in the naive map as well as a "pure" sky 
map produced by ROMAGAL. 

In Figure[l3]a comparison between the naive map and the RO- 
MAGAL map of / = 30° field at 70yt/m is shown. The GLS code is 
capable to remove the 1 // residuals without loosing any signal both 
on bright sources and on diffuse emission. 

In particular, we choose three main proxies: 

• the difference between naive and igls should show only a pat- 
tern due to the l/f noise residuals in the binned map. The pattern 
of this low-frequency noise is recognizable as stripes superimposed 
on the sky signal in the naive map. The stripes are the consequence 
of the l/f noise due to the adopted scanning strategy. 

In Figure [14] we show the difference between igls map and naive 
map. The l/f noise is removed in the igls map but not in the naive. 
The residual stripes due to the low-frequency noise are clearly vis- 
ible. 

• the source fluxes should be the same between the igls and 
naive maps. 

This item is quantified by the map difference, where the pattern 
is due only to the noise without any residual signal, except across 
the very brilliant source. This last effect is discussed on the next 
Section. 

• a statistical analysis of the background noise level should 





PACS / = 


30° field 




Band 


rms igls (MJy/pixel) 


rms naive(MJy/pixel) 


ratio 


70/im 


0.0085 


0.026 


-3.1 


160//m 


0.047 


0.102 


-2.2 




PACS / = 


59° field 




Band 


rms igls (MJy/pixel) 


rms naive(MJy/pixel) 


ratio 


70/zm 


0.004545 


0.02208 


-4.9 


160^m 


0.01899 


0.03586 


- 1.9 



Table 3. rms compared from GLS and naive map on both the SDP observa- 
tions, for PACS bands, measured on a background region of 50x50 pixels. 
In the last column is reported the ratio between the naive and GLS rms. 



show a decrease of the rms value in the igls with respect to the 
naive map. 

In Tables [3] and HI we report the rms residuals of the PACS and 
SPIRE maps respectively, calculated in a diffuse emission area of 
each map. Since the flux between the maps is conserved, a decreas- 
ing of the rms noise level assures an increasing of S/N ratio in the 
ROMAGAL maps. 

The ratio between the naive and igls rms shows an improvement 
of a factor ~ 2 - 5 in the PACS ROMAGAL maps, and a factor 
~ 1 - 2 in the SPIRE case. The difference is mostly due to an 
intrinsically different l/f noise level. 




SPIRE / = 30° field 

Band rms igls (MJy/beam) rms naive(MJy/beam) ratio 

250//m 0.1749 0.2868 ~ 1.6 

350/mi 0.1569 0.2302 -1.5 

500//m 0.2659 0.4065 ~ 1.5 

SPIRE / = 59° field 

Band rms igls (MJy/beam) rms naive(MJy/beam) ratio 

250//m 0.09857 0.1123 - 1.1 

350/im 0.0734 0.08164 ~ 1.1 

500//m 0.1073 0.2101 ~ 1.9 



Table 4. rms compared from GLS and naive map on both the SDP observa- 
tions, for SPIRE bands, measured on a background region of 50x50 pixels. 
In the last column is reported the ratio between the naive and GLS rms. The 
1 // noise is less evident in the SPIRE bolometers with respect to the PACS 
one, but its effect is still remarkable. 



5.3 Consistency of data analysis assumptions 

One of the assumptions of ROMAGAL, as well as of all Fourier 
based GLS map making codes, is that the underlying sky signal is 
constant within a chosen resolution element. If this is not the case, 
artifacts (stripes) will be generat ed in the final map, co ntributing to 
the so called pixelization noise dPoutanen et al]l2006h . In the case 
of Hi-GAL the situation is complicated by several effects: 

• On-board coaddition of samples: each PACS 70yi/m (160 iim) 
frame is the result of an on-board average of eight (four) consecu- 
tive frames reducing the effective sampling frequency of the instru- 
ment (see Section [5Jl ). Thus, sky signal is low-pass filtered by only 
partially effective data- windowing, rather then telescope response, 
leaving room for signal aliasing. 

The map making code is quite sensitive to aliasing since it works 
in Fourier space. The situation is worsened by the large dynamic 
range of the Hi-GAL fields, especially when scanning across bright 
sources. 

• Time bolometer response induces signal distortions along the 
scan. While within HIPE the SPIRE detector response is decon- 
volved from the data (Griffin et al]l2008h . the same is not true for 
PACS. Redundancy in each pixel is obtained by scans coming from 
different directions, thus the effect contributes further signal mis- 
match at the pixel level. 

• P ointing error: as analyzed in detail in iGarcia Lario et al.l 
d2007h the pointing performance of Herschel, which mean the ca- 
pability of assign coordinates to each sample in a given refer- 
ence frame, can be affected by several pointing error effects; the 
main contributor is due to the angular separation between the de- 
sired direction and the actual instantaneous direction, based on the 
position-dependent bias within the star- trackers. 

The Herschel AOCS goal is that the mismatch between real co- 
ordinates and assign ed coordinates along the scan-leg is smaller 
than 1.2" at 1 sigma jGarcia Lario et al.ll2007l) . So that a 2cr event 
becomes significant compared to the PSF of the PACS blue band. 

All of the above effects challenge the basic assumptions (no 
signal aliasing and no sub-pixel effect) under which ROMAGAL 
works. Our simulations suggest that signal aliasing contribute sig- 



Figure 15. Zoom around a compact sources for the PACS blue optimal map. 
The striping dragged along the scan directions are remarkable. 

nificantly more than the other two. The net result on maps is the 
striping of bright sources in the GLS maps. An example is shown 
in Figure Q3] for the PACS blue band. 

It is important to notice that the stripes are present only around 
the point-like sources (where - of course - the signal aliasing is 
more evident), regardless of their magnitude. However, magnitude 
influences the amplitude of the stripes. For a source peak which is 
only 10 times higher that the rms background, the intensity of the 
stripes are within lcr from the background dispersion. For the more 
intense sources, the stripes magnitude in the pixels surrounding the 
PSF shape can be 100cr times away from the background value. 

Since these stripes are produced within the GLS solver, which 
performs repetitive convolutions along the scan directions, but do 
not affect the naive map the obvious workaround is to implement 
dedicated map making which considers a different analysis around 
the bright sources. 

However, the detailed accounting of the above effects and the 
enhanced map making strategies to a ddress them will be the subject 
of a forthcoming paper lPiazzo et al.ld201lh . 



6 SUMMARY AND CONCLUSIONS 

This paper describes in detail all steps of the processing of Herschel 
data, from originally downloaded frames to high-quality maps, 
used in the framework of Hi-GAL survey. Hi-GAL data are taken in 
fast scan mode (6(y7sec ) and simultaneously by PACS and SPIRE 
(Parallel mode). We test our pipeline reducing data from the Sci- 
ence Demonstration Phase and present results, taking as proxy for 
the quality of the final images their comparison with naive maps. 

We divided data processing into two distinct phases: prepro- 
cessing and mapmaking. Pre-processing aims to accurately remove 
systematics and random effects in the data, in order to prepare them 
for the ROMAGAL map making algorithm, which implements the 
minimum variance GLS approach in order to minimize the noise 
into the data. It turns out that NCR is a fundamental step in the 
pre-processing because ROMAGAL, as an FFT mapmaking code 
needs continuous and regular data time series as input. 

Noise residuals in the diffuse emission of the two test fields 
(SDP Hi-GAL data, two 2° x 2° tiles centered on the Galactic plane 
at 1 = 30° and 1 = 59°) show that we obtain optimal maps, getting 
rid of glitches and systematic drifts, as well as minimizing the 1/ 
f and white noise component. The remaining effects, which do not 
affect the overall quality of the maps except across the bright source 
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on the PACS 70//m maps, are under investigation that will appear 
in a dedicated publication to be available shortly. 
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